function [out, vars, vectors] = GPOPS_plot_results(output,setup)
% computes variables based on GPOPS result, including various
% statistics, and terms from the optimal tax formula. Also does
% some plotting, but that is not the main purpose. Stores
% transformed results.
    
tic
%% assign variables
    auxdata = setup.auxdata;
    
    %% Parameters
    inpscale = auxdata.settings.inpscale;
    out.L_M = inpscale.L_M.* output.result.solution.parameter(1,1);
    out.L_R = inpscale.L_R.* output.result.solution.parameter(1,2);
    out.L_C = inpscale.L_C.* output.result.solution.parameter(1,3);
    out.K_B_M = inpscale.K_B_M.* output.result.solution.parameter(1,4);
    out.K_B_R = inpscale.K_B_R.* output.result.solution.parameter(1,5);
    out.K_B_C = inpscale.K_B_C.* output.result.solution.parameter(1,6);
    out.K_E = inpscale.K_E.* output.result.solution.parameter(1,7);
    out.K_S = inpscale.K_S.* output.result.solution.parameter(1,8);
    out.b = inpscale.transfer.* output.result.solution.parameter(1,9);

    out.K_B = out.K_B_M + out.K_B_R + out.K_B_C;
    
    % derived factor prices
    gpops_fun = auxdata.gpops_fun;

    out.Y = gpops_fun.prod_fun.Y(out);
    out.Y_M = gpops_fun.prod_fun.Y_M(out);
    out.Y_R = gpops_fun.prod_fun.Y_R(out);
    out.Y_C = gpops_fun.prod_fun.Y_C(out);
    out.Y_B_M = gpops_fun.prod_fun.Y_B_M(out);
    out.Y_B_R = gpops_fun.prod_fun.Y_B_R(out);
    out.Y_B_C = gpops_fun.prod_fun.Y_B_C(out);
    out.Y_E = gpops_fun.prod_fun.Y_E(out);
    out.Y_S = gpops_fun.prod_fun.Y_S(out); 

    
    out.wtilde_orig = output.result.interpsolution.phase.time;
    out.wtilde_noninterp = output.result.solution.phase.time;
    
    % interpolate
    interp_points = 10000; % choosing many points for interpolation
                           % of objects (not in the programming
                           % sense) which are used to compute the
                           % terms in the optimal tax formula. It
                           % seems that the error in these terms
                           % goes down with a large number of
                           % points. Interpolation takes a long
                           % time, though.
    
    out.wtilde = linspace(out.wtilde_orig(1),out.wtilde_orig(end),interp_points)';
    
    if auxdata.settings.log_w == 1
        if auxdata.settings.standardize == 1
            out.w_orig = exp(out.wtilde_orig.*auxdata.sigma_fit+auxdata.mu_fit);
            out.w = exp(out.wtilde.*auxdata.sigma_fit+auxdata.mu_fit);
            out.w_noninterp = exp(out.wtilde_noninterp.*auxdata.sigma_fit+auxdata.mu_fit);
        else
            out.w_orig = exp(out.wtilde_orig);
            out.w = exp(out.wtilde);
            out.w_noninterp = exp(out.wtilde_noninterp);
        end
    end
    
    out.w0 = auxdata.w0;
    out.wf = auxdata.wf;

    out.u = output.result.interpsolution.phase.state .* inpscale.u;
    out.u_noninterp = output.result.solution.phase.state .* inpscale.u;
    out.l = output.result.interpsolution.phase.control(:,1) .* inpscale.l;
    out.l_noninterp = output.result.solution.phase.control(:,1) .*inpscale.l;

    % multipliers
    out.mu = -output.result.solution.phase.integralmultipliers(5);
    out.xi_M = -output.result.solution.phase.integralmultipliers(6)./out.mu*out.Y_M;
    out.xi_R = -output.result.solution.phase.integralmultipliers(7)./out.mu*out.Y_R;
    out.xi_C = -output.result.solution.phase.integralmultipliers(8)./out.mu*out.Y_C;
    out.eta = output.result.interpsolution.phase.costate/out.mu;
    
    %% finer interpolation
    out.l_orig = out.l;
    out.u_orig = out.u;
    out.eta_orig = out.eta;
    
    out.l = interp1(out.w_orig,out.l_orig,out.w);
    out.u = interp1(out.w_orig,out.u_orig,out.w);
    out.eta = interp1(out.w_orig,out.eta_orig,out.w);
    
    %% assign parameters
    epsilon = gpops_fun.util_fun.epsilon;
    k = 1/epsilon;
    
    xi = gpops_fun.util_fun.xi;
    %% construct c
    out.c = 1./xi.*out.l.^(1+1/epsilon)/(1+1/epsilon) + out.u;
    
    %% compute derived variables
    out.mtax = 1-1./out.w.*1./xi.*out.l.^k;

    out.y = out.l.*out.w;
    out.tax = out.y - out.c;
    out.avg_tax = out.tax./out.y;
    
    % total mass
    out.total_mass = auxdata.gpops_fun.skill_dist.mass;
    
    % capital income and capital expenditure
    out.capital_income = out.K_B_M .* out.Y_B_M + ...
                out.K_B_R .* out.Y_B_R + ...
                out.K_B_C .* out.Y_B_C + ...
                out.K_E .* out.Y_E + ...
                out.K_S .* out.Y_S;

    out.q_B = gpops_fun.q_B;
    out.q_E = gpops_fun.q_E;
    out.q_S = gpops_fun.q_S;
    out.drop = auxdata.drop;
    
    out.capital_expenditure = out.K_B_M .*out.q_B + ...
                out.K_B_R .* out.q_B + ...
                out.K_B_C .* out.q_B + ...
                out.K_E .* out.q_E + ...
                out.K_S .* out.q_S;
    
    % tax revenue and value of resource constraints
    out.labor_income = out.L_M * out.Y_M + ...
        out.L_R * out.Y_R + ...
        out.L_C * out.Y_C;

    % income shares
    out.labor_income_share = out.labor_income/out.Y;
    out.robot_income = out.K_B * out.Y_B_M; % requires equal
                                            % returns to all robot capital
    out.robot_income_share = out.robot_income/out.Y;
    out.labor_share_in_labor_plus_robot_income = out.labor_income/ ...
        (out.labor_income + out.robot_income);
    
    % bounds check
    fprintf('solution - lower bound, should be positive\n')
    check_lb = output.result.solution.parameter - setup.bounds.parameter.lower
    
    fprintf('upper bound - solution, should be positive\n')
    check_ub = setup.bounds.parameter.upper - output.result.solution.parameter
    
    % taxes
    out.tau_B = out.Y_B_M/out.q_B - 1;
    out.tau_E = out.Y_E/out.q_E - 1;
    out.tau_S = out.Y_S/out.q_S - 1;
    
    out.tau_return_B = 1 - 1/(1+out.tau_B);
    out.tau_return_E = 1 - 1/(1+out.tau_E);
    out.tau_return_S = 1 - 1/(1+out.tau_S);

    out.fval = output.result.objective * inpscale.welfare;
    out.welfare =  -output.result.objective * inpscale.welfare;
    
    out.welfare_M = output.result.solution.phase.integral(1);
    out.welfare_R = output.result.solution.phase.integral(2);
    out.welfare_C = output.result.solution.phase.integral(3);
    out.welfare_nonpart = output.result.solution.phase.integral(4);
    
    %% population shares, participation and employment
    w_lb = out.w(1);
    w_ub = out.w(end);
    
    % approximate phi_tilde
    out.phi_tilde_w = out.u - out.b;
    
    phi_tilde_fun = @(w) interp1(out.w,out.phi_tilde_w,w,'spline','extrap');
    
    w_range = w_lb:w_ub;
    figure(7)
    plot(out.w,out.phi_tilde_w,w_range,phi_tilde_fun(w_range),'--')
    title('phi tilde')
    legend('actual','approximated')
    
    
    f_M = @(w) gpops_fun.skill_dist.f_M(w,out.Y_M,out.Y_R,out.Y_C);
    f_R = @(w) gpops_fun.skill_dist.f_R(w,out.Y_M,out.Y_R,out.Y_C);
    f_C = @(w) gpops_fun.skill_dist.f_C(w,out.Y_M,out.Y_R,out.Y_C);

    out.f_M = f_M(out.w);
    out.f_R = f_R(out.w);
    out.f_C = f_C(out.w);
    out.f = out.f_M + out.f_R + out.f_C;
    
    out.pop_M = integral(f_M,w_lb,w_ub);
    out.pop_R = integral(f_R,w_lb,w_ub);
    out.pop_C = integral(f_C,w_lb,w_ub);
    
    if auxdata.settings.fix_participation
        
        pi_M = auxdata.settings.pi_M;
        pi_R = auxdata.settings.pi_R;
        pi_C = auxdata.settings.pi_C;
    
        h_M = @(w) gpops_fun.skill_dist.f_M(w,out.Y_M,out.Y_R,out.Y_C).*pi_M;
        h_R = @(w) gpops_fun.skill_dist.f_R(w,out.Y_M,out.Y_R,out.Y_C).*pi_R;
        h_C = @(w) gpops_fun.skill_dist.f_C(w,out.Y_M,out.Y_R,out.Y_C).*pi_C;
    
    else
        h_M = @(w) gpops_fun.skill_dist.f_M(w,out.Y_M,out.Y_R,out.Y_C).*gpops_fun.part_dist.G_M(phi_tilde_fun(w));
        h_R = @(w) gpops_fun.skill_dist.f_R(w,out.Y_M,out.Y_R,out.Y_C).*gpops_fun.part_dist.G_R(phi_tilde_fun(w));
        h_C = @(w) gpops_fun.skill_dist.f_C(w,out.Y_M,out.Y_R,out.Y_C).*gpops_fun.part_dist.G_C(phi_tilde_fun(w));

    end
    
    out.h_M = h_M(out.w);
    out.h_R = h_R(out.w);
    out.h_C = h_C(out.w);
    out.h = out.h_M + out.h_R + out.h_C;    

    out.F_M = zeros(size(out.w));
    out.F_R = zeros(size(out.w));
    out.F_C = zeros(size(out.w));
    out.H_M = zeros(size(out.w));
    out.H_R = zeros(size(out.w));
    out.H_C = zeros(size(out.w));

    % can use G_M here since participation costs are participation-independent
    out.pi_w = gpops_fun.part_dist.G_M(phi_tilde_fun(out.w)); 
    
    for i=1:length(out.H_M)
        out.F_M(i) = integral(f_M,w_lb,out.w(i));
        out.F_R(i) = integral(f_R,w_lb,out.w(i));
        out.F_C(i) = integral(f_C,w_lb,out.w(i));
        
        out.H_M(i) = integral(h_M,w_lb,out.w(i));
        out.H_R(i) = integral(h_R,w_lb,out.w(i));
        out.H_C(i) = integral(h_C,w_lb,out.w(i));
    end
    
    out.F = out.F_M + out.F_R + out.F_C;
    out.H = out.H_M + out.H_R + out.H_C;
    
    out.part_M = integral(h_M,w_lb,w_ub);
    out.part_R = integral(h_R,w_lb,w_ub);
    out.part_C = integral(h_C,w_lb,w_ub);
    
    out.part = out.part_M + out.part_R + out.part_C;
    
    % average participation rates by occ.
    out.pi_M = out.part_M/out.pop_M;
    out.pi_R = out.part_R/out.pop_R;
    out.pi_C = out.part_C/out.pop_C;
    
    out.h_M_cond_part = out.h_M / out.part_M;
    out.h_R_cond_part = out.h_R / out.part_R;
    out.h_C_cond_part = out.h_C / out.part_C;
    out.h_cond_part = out.h/out.part;

    out.H_M_cond_part = out.H_M / out.part_M;
    out.H_R_cond_part = out.H_R / out.part_R;
    out.H_C_cond_part = out.H_C / out.part_C;
    out.H_cond_part = out.H/out.part;
    
    out.emp_M = out.part_M/out.part;
    out.emp_R = out.part_R/out.part;
    out.emp_C = out.part_C/out.part;
    
    out.average_wage_M = integral(@(w) w.*h_M(w),w_lb,w_ub)*1/out.emp_M*1/out.part;
    out.average_wage_R = integral(@(w) w.*h_R(w),w_lb,w_ub)*1/out.emp_R*1/out.part;
    out.average_wage_C = integral(@(w) w.*h_C(w),w_lb,w_ub)*1/out.emp_C*1/out.part;
    
    out.total_cons_M = output.result.solution.phase.integral(9);
    out.total_cons_R = output.result.solution.phase.integral(10);
    out.total_cons_C = output.result.solution.phase.integral(11);
    
    % average income by occupation
    out.avg_inc_M = out.L_M*out.Y_M/out.emp_M;
    out.avg_inc_R = out.L_R*out.Y_R/out.emp_R;
    out.avg_inc_C = out.L_C*out.Y_C/out.emp_C;
    
    %% formula terms
    % here I compute the various terms in the optimal tax
    % formulas. Many of the terms are only here for diagnostic
    % reasons and are not plotted in the paper.
    
    % parts of IC effects
    dV = diff(out.u);
    dw = diff(out.w);
    dV_dw = dV./dw;
    out.V_prime = [dV_dw(1);dV_dw];
    out.dw = [dw(1);dw];

    dh_M_h = diff(out.h_M./out.h);
    dh_R_h = diff(out.h_R./out.h);
    dh_C_h = diff(out.h_C./out.h);
    dh_M_h_dw = dh_M_h./dw;
    dh_R_h_dw = dh_R_h./dw;
    dh_C_h_dw = dh_C_h./dw;

    out.dh_M_h_dw = [dh_M_h_dw(1);dh_M_h_dw];
    out.dh_R_h_dw = [dh_R_h_dw(1);dh_R_h_dw];
    out.dh_C_h_dw = [dh_C_h_dw(1);dh_C_h_dw];

    % IC effects
    out.I_M = sum(out.eta .* out.V_prime .* out.w .* out.dh_M_h_dw .* out.dw);
    out.I_R = sum(out.eta .* out.V_prime .* out.w .* out.dh_R_h_dw .* out.dw);
    out.I_C = sum(out.eta .* out.V_prime .* out.w .* out.dh_C_h_dw .* out.dw);

    % parts of C effects
    dl = diff(out.l);
    dl_dw = dl./dw;
    out.l_prime = [dl_dw(1);dl_dw];

    out.C_RM = 1/out.Y_M * sum(out.l_prime .* out.w.^2 .* out.h_R .* ...
                               out.h_M ./ out.h .* out.dw);

    out.C_CM = 1/out.Y_M * sum(out.l_prime .* out.w.^2 .* out.h_C .* ...
                               out.h_M ./ out.h .* out.dw);

    out.C_MR = 1/out.Y_R * sum(out.l_prime .* out.w.^2 .* out.h_M .* ...
                               out.h_R ./ out.h .* out.dw);

    out.C_CR = 1/out.Y_R * sum(out.l_prime .* out.w.^2 .* out.h_C .* ...
                               out.h_R ./ out.h .* out.dw);

    out.C_MC = 1/out.Y_C * sum(out.l_prime .* out.w.^2 .* out.h_M .* ...
                               out.h_C ./ out.h .* out.dw);

    out.C_RC = 1/out.Y_C * sum(out.l_prime .* out.w.^2 .* out.h_R .* ...
                               out.h_C ./ out.h .* out.dw);

    out.C_MM = - out.C_MR - out.C_MC;
    out.C_RR = - out.C_RM - out.C_RC;
    out.C_CC = - out.C_CM - out.C_CR;

    out.xi_C_iM_sum = out.xi_M * out.C_MM + out.xi_R * out.C_RM + out.xi_C * out.C_CM;
    out.xi_C_iC_sum = out.xi_M * out.C_MC + out.xi_R * out.C_RC + out.xi_C * out.C_CC;
    
    out = compute_elasticities(out, 0.01, setup);

    out.LHS = out.tau_B * out.K_B * out.q_B;
    out.RHS_no_part = out.eps_Y_M_Y_R_K_B * (out.I_M + ...
                                               out.xi_M * out.C_MM  + ...
                                               out.xi_R * out.C_RM  + ...
                                               out.xi_C * out.C_CM) + ...
        out.eps_Y_C_Y_R_K_B * (out.I_C + ...
                                 out.xi_M * out.C_MC  + ...
                                 out.xi_R * out.C_RC  + ...
                                 out.xi_C * out.C_CC);

    out.part_residual = out.LHS - out.RHS_no_part;

    %% terms for participation effects    
    dpi_w = diff(out.pi_w);
    dpi_w_dw = dpi_w./dw;
    out.pi_prime_w = [dpi_w_dw(1); dpi_w_dw];
    
    dc_dw = diff(out.c)./dw;
    out.c_prime = [dc_dw(1); dc_dw];
    
    out.Pi_M = cumsum(out.w .* out.pi_prime_w .* out.f_M .* out.dw);
    out.Pi_R = cumsum(out.w .* out.pi_prime_w .* out.f_R .* out.dw);
    out.Pi_C = cumsum(out.w .* out.pi_prime_w .* out.f_C .* out.dw);
    
    dPi_M_wh_dw = diff(out.Pi_M./(out.w .* out.h))./dw;
    dPi_R_wh_dw = diff(out.Pi_R./(out.w .* out.h))./dw;
    dPi_C_wh_dw = diff(out.Pi_C./(out.w .* out.h))./dw;
    
    out.dPi_M_wh_dw = [dPi_M_wh_dw(1); dPi_M_wh_dw];
    out.dPi_R_wh_dw = [dPi_R_wh_dw(1); dPi_R_wh_dw];
    out.dPi_C_wh_dw = [dPi_C_wh_dw(1); dPi_C_wh_dw];
    
    out.I_Pi_M = sum(out.eta .* out.V_prime .* out.w .* out.dPi_M_wh_dw .* out.dw);
    out.I_Pi_R = sum(out.eta .* out.V_prime .* out.w .* out.dPi_R_wh_dw .* out.dw);
    out.I_Pi_C = sum(out.eta .* out.V_prime .* out.w .* out.dPi_C_wh_dw .* out.dw);
    
    out.F_Pi_M = sum( (out.l .* out.w - out.c + out.b)  ...
                 .* out.pi_prime_w./out.pi_w .* out.Pi_M .* out.dw);

    out.F_Pi_R = sum( (out.l .* out.w - out.c + out.b)  ...
                 .* out.pi_prime_w./out.pi_w .* out.Pi_R .* out.dw);
    
    out.F_Pi_C = sum( (out.l .* out.w - out.c + out.b)  ...
                 .* out.pi_prime_w./out.pi_w .* out.Pi_C .* out.dw);
    
    out.B_Pi_M = sum((out.l_prime .* out.w - out.c_prime) .* out.Pi_M ...
                 .* out.dw);

    out.B_Pi_R = sum((out.l_prime .* out.w - out.c_prime) .* out.Pi_R ...
                 .* out.dw);
    
    out.B_Pi_C = sum((out.l_prime .* out.w - out.c_prime) .* out.Pi_C ...
                 .* out.dw);    
    
    C_Pi_B_sum = out.dY_M_dK_B/out.Y_M .* out.Pi_M + ...
                 out.dY_R_dK_B/out.Y_R .* out.Pi_R + ...
                 out.dY_C_dK_B/out.Y_C .* out.Pi_C;
    
    C_Pi_M = out.xi_M * 1/out.Y_M * sum(out.l_prime .* out.w .* ...
                                        out.h_M./out.h .* C_Pi_B_sum ...
                                        .* out.dw);

    C_Pi_R = out.xi_R * 1/out.Y_R * sum(out.l_prime .* out.w .* ...
                                        out.h_R./out.h .* C_Pi_B_sum ...
                                        .* out.dw);    
    
    C_Pi_C = out.xi_C * 1/out.Y_C * sum(out.l_prime .* out.w .* ...
                                        out.h_C./out.h .* C_Pi_B_sum ...
                                        .* out.dw);    
    
    out.C_Pi_B = out.K_B * (C_Pi_M + C_Pi_R + C_Pi_C);
    
    % terms for participation shift
    out.dF_M_dK_B = - out.dY_M_dK_B/out.Y_M .* out.w .* out.f_M;
    out.dF_R_dK_B = - out.dY_R_dK_B/out.Y_R .* out.w .* out.f_R;
    out.dF_C_dK_B = - out.dY_C_dK_B/out.Y_C .* out.w .* out.f_C;
    out.dF_dK_B = out.dF_M_dK_B + out.dF_R_dK_B + out.dF_C_dK_B;
    
    
    out.dH_M_dK_B = out.pi_w .* out.dF_M_dK_B + out.dY_M_dK_B/out.Y_M .* out.Pi_M;
    out.dH_R_dK_B = out.pi_w .* out.dF_R_dK_B + out.dY_R_dK_B/out.Y_R .* out.Pi_R;
    out.dH_C_dK_B = out.pi_w .* out.dF_C_dK_B + out.dY_C_dK_B/out.Y_C .* out.Pi_C;
    
    out.dH_dK_B = out.dH_M_dK_B + out.dH_R_dK_B + out.dH_C_dK_B;
    
    
    out.S_Pi_B = out.K_B * ( out.xi_M * 1/out.Y_M * sum(out.l .* out.w .* out.pi_prime_w .* (out.f_M./out.h .* out.dH_dK_B - out.dF_M_dK_B) .* out.dw) + ...
                             out.xi_R * 1/out.Y_R * sum(out.l .* out.w .* out.pi_prime_w .* (out.f_R./out.h .* out.dH_dK_B - out.dF_R_dK_B) .* out.dw) + ...
                             out.xi_C * 1/out.Y_C * sum(out.l .* out.w .* out.pi_prime_w .* (out.f_C./out.h .* out.dH_dK_B - out.dF_C_dK_B) .* out.dw));
    
    out.W_Pi_M = 1/out.mu*sum(welfare_shift_M(gpops_fun, out) .* out.dw);
    out.W_Pi_R = 1/out.mu*sum(welfare_shift_R(gpops_fun, out) .* out.dw);
    out.W_Pi_C = 1/out.mu*sum(welfare_shift_C(gpops_fun, out) .* out.dw);
    
    out.part_effect_term_1_M = out.eps_Y_M_K_B * (out.W_Pi_M + ...
                                                    out.I_Pi_M + ...
                                                    out.F_Pi_M + out.B_Pi_M );

    out.part_effect_term_1_R = out.eps_Y_R_K_B * (out.W_Pi_R + ...
                                                    out.I_Pi_R + ...
                                                    out.F_Pi_R + ...
                                                    out.B_Pi_R );
    
    out.part_effect_term_1_C = out.eps_Y_C_K_B * (out.W_Pi_C + ...
                                                    out.I_Pi_C + ...
                                                    out.F_Pi_C + ...
                                                    out.B_Pi_C );
    
    out.part_effect = - out.part_effect_term_1_M ...
                      - out.part_effect_term_1_R ...
                      - out.part_effect_term_1_C ...
                      - out.C_Pi_B ...
                      - out.S_Pi_B;
    
    out.I_Pi_parts = out.eps_Y_M_K_B*out.I_Pi_M + ...
                     out.eps_Y_R_K_B*out.I_Pi_R + ...
                     out.eps_Y_C_K_B*out.I_Pi_C;
    
    out.W_Pi_parts = out.eps_Y_M_K_B*out.W_Pi_M + ...
                     out.eps_Y_R_K_B*out.W_Pi_R + ...
                     out.eps_Y_C_K_B*out.W_Pi_C;
    
    out.F_Pi_parts = out.eps_Y_M_K_B*out.F_Pi_M + ...
                     out.eps_Y_R_K_B*out.F_Pi_R + ...
                     out.eps_Y_C_K_B*out.F_Pi_C;
    
    out.B_Pi_parts = out.eps_Y_M_K_B*out.B_Pi_M + ...
                     out.eps_Y_R_K_B*out.B_Pi_R + ...
                     out.eps_Y_C_K_B*out.B_Pi_C;    

    
    out.RHS_with_effects = out.RHS_no_part + out.part_effect;
    
    %% FOC in early step
    out.mrs = 1/xi*out.l.^(1./epsilon);
    
    out.foc_row_1 = -sum(welfare_term_foc(gpops_fun, out) .* out.dw);
    out.foc_row_2 = -out.mu * ( out.xi_M * out.dY_M_dK_B/out.Y_M*out.L_M ...
                              + out.xi_R * out.dY_R_dK_B/out.Y_R*out.L_R ...
                              + out.xi_C * out.dY_C_dK_B/out.Y_C* ...
                                out.L_C);
    out.foc_row_3 = -out.mu * (   out.xi_M * 1/out.Y_M * sum( (out.l_prime.*out.w + out.l) .* out.dH_M_dK_B .* out.dw)...
                                + out.xi_R * 1/out.Y_R * sum( (out.l_prime.*out.w + out.l) .* out.dH_R_dK_B .* out.dw)...
                                + out.xi_C * 1/out.Y_C * sum( (out.l_prime.*out.w + out.l) .* out.dH_C_dK_B .* out.dw));
    
   
    
    out.foc_row_4 = -out.mu .* sum( (out.l_prime .* out.w + out.l - ...
                                     out.V_prime - out.mrs.*out.l_prime) .* out.dH_dK_B .* out.dw);
    
    out.foc_row_5 = -out.mu * out.b * sum(out.pi_prime_w .* out.dF_dK_B ...
                                          .* out.dw);
    
    out.foc_row_6 = out.mu * ( out.Y_B_M - out.q_B + ...
                                out.K_B * out.dY_B_dK_B ...
                              + out.K_E * out.dY_E_dK_B ...
                              + out.K_S * out.dY_S_dK_B);
    
    
    
    % B_1
    B_1_row_1_tmp = welfare_term_B_1(gpops_fun, out);
    B_1_row_2_tmp = - Psi(gpops_fun,out.b).*out.h;
    B_1_row_3_tmp = out.mu.*(  out.xi_M/out.Y_M.*out.l.*out.w.*out.h_M ...
                             + out.xi_R/out.Y_R.*out.l.*out.w.*out.h_R ...
                             + out.xi_C/out.Y_C.*out.l.*out.w.*out.h_C);
    B_1_row_4_tmp = out.mu.*(out.l.*out.w - out.c + out.b).*out.h;
    
    % boundary wage derivatives, assuming that top wage is earned
    % in C and bottom wage in M
    dw_lb_dK_B = out.w(1)*out.dY_M_dK_B/out.Y_M;
    dw_ub_dK_B = out.w(end)*out.dY_C_dK_B/out.Y_C;
    
    B_1_row_1 = dw_ub_dK_B * B_1_row_1_tmp(end) - dw_lb_dK_B * B_1_row_1_tmp(1);
    B_1_row_2 = dw_ub_dK_B * B_1_row_2_tmp(end) - dw_lb_dK_B * B_1_row_2_tmp(1);
    B_1_row_3 = dw_ub_dK_B * B_1_row_3_tmp(end) - dw_lb_dK_B * B_1_row_3_tmp(1);
    B_1_row_4 = dw_ub_dK_B * B_1_row_4_tmp(end) - dw_lb_dK_B * B_1_row_4_tmp(1);
    
    out.B_1 = B_1_row_1 + B_1_row_2 + B_1_row_3 + B_1_row_4;
    
    % B_2
    tmp = welfare_term_B_2(gpops_fun, out);
    out.B_2_row_1 = tmp(end) - tmp(1);
    tmp = out.mu * (  out.xi_M / out.Y_M .* out.l .* out.w .* out.dH_M_dK_B ...
                    + out.xi_R / out.Y_R .* out.l .* out.w .* out.dH_R_dK_B ...
                    + out.xi_C / out.Y_C .* out.l .* out.w .* out.dH_C_dK_B);
    out.B_2_row_2 = tmp(end) - tmp(1);
    out.B_2_row_3 = out.mu * ((out.l(end) * out.w(end) - out.c(end))*out.dH_dK_B(end)...
                             -(out.l(1) * out.w(1) - out.c(1))* ...
                              out.dH_dK_B(1));
    tmp = (out.mu * out.b - Psi(gpops_fun,out.b)) .* out.pi_w .* out.dF_dK_B;
    out.B_2_row_4 = tmp(end) - tmp(1);
    
    out.B_2 = out.B_2_row_1 + out.B_2_row_2 + out.B_2_row_3 + ...
              out.B_2_row_4;
    

    
    % putting everything together
    out.foc = out.foc_row_1 + out.foc_row_2 + out.foc_row_3 + ...
              out.foc_row_4 + out.foc_row_5 + out.foc_row_6 + out.B_1 ...
              + out.B_2;
    
    out.foc_RHS = - out.foc / out.mu * out.K_B + out.tau_B*out.K_B*out.q_B; 
    out.robot_tax_rev = out.tau_B*out.K_B*out.q_B;
    
    out.robot_tax_base = out.K_B*out.q_B;
    %% equation (231)
    
    out.eq231_row_1 = out.mu .* sum((out.l.*out.w - out.c + ...
                                     out.b).*out.pi_prime_w./out.pi_w.*out.dH_dK_B.*out.dw);
    out.eq231_row_2 = -out.mu * (  out.dY_M_dK_B/out.Y_M*sum(out.l.*out.Pi_M.*out.dw)...
                                 + out.dY_R_dK_B/out.Y_R*sum(out.l.*out.Pi_R.*out.dw)...
                                 + out.dY_C_dK_B/out.Y_C* ...
                                   sum(out.l.*out.Pi_C.*out.dw)  );
    
    tmp = out.mu .* (out.l.*out.w - out.c + out.b) .* (  out.dY_M_dK_B./out.Y_M.*out.Pi_M...
                                                   + out.dY_R_dK_B./out.Y_R.*out.Pi_R... 
                                                   + ...
                                                     out.dY_C_dK_B./out.Y_C.*out.Pi_C);
    out.eq231_row_3 = tmp(end) - tmp(1);
       
    out.eq231 = out.eq231_row_1 + out.eq231_row_2 + ...
        out.eq231_row_3;
    
    %% equation (237)
    
    out.eq237_row_1 = out.mu .* sum( (out.l.*out.w - out.c + out.b ).*out.pi_prime_w./out.pi_w.*...
                  (  out.dY_M_dK_B/out.Y_M.*out.Pi_M...
                   + out.dY_R_dK_B/out.Y_R.*out.Pi_R...
                   + out.dY_C_dK_B/out.Y_C.*out.Pi_C ).*out.dw);
    
    out.eq237_row_2 = out.mu * sum( (out.l_prime.*out.w - out.c_prime).*...
                  (  out.dY_M_dK_B/out.Y_M.*out.Pi_M...
                   + out.dY_R_dK_B/out.Y_R.*out.Pi_R...
                   + out.dY_C_dK_B/out.Y_C.*out.Pi_C ).*out.dw);
    
    out.eq237 = out.eq237_row_1 + out.eq237_row_2;
    
    %% equation (240)mn
    out.eq240 = out.mu * (  out.dY_M_dK_B/out.Y_M * (out.F_Pi_M + out.B_Pi_M)...
                          + out.dY_R_dK_B/out.Y_R * (out.F_Pi_R + out.B_Pi_R)...
                          + out.dY_C_dK_B/out.Y_C * (out.F_Pi_C + out.B_Pi_C));
    
    
    %% equation (187)
    
    out.OB = sum(welfare_term_OB(gpops_fun, out).*...
             out.V_prime./out.pi_w .* ...
             (  out.dY_M_dK_B/out.Y_M .* out.Pi_M ...
              + out.dY_R_dK_B/out.Y_R .* out.Pi_R ...
              + out.dY_C_dK_B/out.Y_C .* out.Pi_C) .* out.dw);
    
    deta_dw = diff(out.eta)./dw;
    out.eta_prime = [deta_dw(1);deta_dw];
    
    dMRS_l_dw = diff(out.mrs.*out.l)./dw;
    out.dMRS_l_dw = [dMRS_l_dw(1);dMRS_l_dw];
    
    out.IC = -out.mu * sum( (out.eta .* out.dMRS_l_dw + out.w .* out.V_prime .* out.eta_prime) .* ...
                            out.dH_dK_B./(out.w .* out.h) .* out.dw);
    
    out.LC_row_1 = out.mu .* sum(out.l .* out.w .* out.pi_prime_w.*...
                                 (  out.xi_M/out.Y_M .* out.f_M./out.h.*out.dH_dK_B...
                                  + out.xi_R/out.Y_R .* out.f_R./out.h.*out.dH_dK_B...
                                  + out.xi_C/out.Y_C .* out.f_C./out.h.*out.dH_dK_B ).*out.dw);
   
    
    out.LC_row_2 = out.mu .* sum(out.l_prime .* out.w .* out.dH_dK_B.*...
                                 (  out.xi_M/out.Y_M .* out.h_M./out.h...
                                  + out.xi_R/out.Y_R .* out.h_R./out.h...
                                  + out.xi_C/out.Y_C .* out.h_C./out.h).*out.dw);
    
    tmp = out.mu.*out.l.*out.w.*(  out.xi_M/out.Y_M.*out.dY_M_dK_B/out.Y_M.*out.Pi_M...
                                 + out.xi_R/out.Y_R.*out.dY_R_dK_B/out.Y_R.*out.Pi_R...
                                 + out.xi_C/out.Y_C.*out.dY_C_dK_B/out.Y_C.*out.Pi_C);
    
    out.LC_row_3 = tmp(end);
    
    out.LC_row_4 = -out.mu .* (  out.xi_M/out.Y_M * (sum( (out.l_prime.*out.w+out.l).*out.dH_M_dK_B.*out.dw) + out.dY_M_dK_B*out.L_M)...
                               + out.xi_R/out.Y_R * (sum( (out.l_prime.*out.w+out.l).*out.dH_R_dK_B.*out.dw) + out.dY_R_dK_B*out.L_R)... 
                               + out.xi_C/out.Y_C * (sum( (out.l_prime.*out.w+out.l).*out.dH_C_dK_B.*out.dw) + out.dY_C_dK_B*out.L_C)); 

    out.LC = out.LC_row_1 + out.LC_row_2 + out.LC_row_3 + out.LC_row_4;

    out.RC_row_1 = out.mu * sum( (out.l.*out.w - out.c + out.b).*out.pi_prime_w./out.pi_w.*out.dH_dK_B.*out.dw );
    
    out.RC_row_2 = -out.mu * sum( out.l.*(  out.dY_M_dK_B/out.Y_M.*out.Pi_M...
                                          + out.dY_R_dK_B/out.Y_R.*out.Pi_R...
                                          + out.dY_C_dK_B/out.Y_C.*out.Pi_C).*out.dw );
    
    out.RC_row_3 = -out.mu * out.b * sum(out.pi_prime_w.*out.dF_dK_B.* out.dw);
    
    tmp = out.mu * (out.l.*out.w - out.c) .* ...
          (  out.dY_M_dK_B/out.Y_M.*out.Pi_M...
           + out.dY_R_dK_B/out.Y_R.*out.Pi_R...
           + out.dY_C_dK_B/out.Y_C.*out.Pi_C );
    
    out.RC_row_4 = tmp(end);
    
    out.RC = out.RC_row_1 + out.RC_row_2 + out.RC_row_3 + out.RC_row_4;
    
    out.eq187 = out.OB + out.IC + out.LC + out.RC + out.mu*(out.Y_B_M ...
                                                      - out.q_B);
    
    %% different parts of collected terms
    
    % OB in eq (194)
    out.eq194 = out.mu * (  out.dY_M_dK_B/out.Y_M * out.W_Pi_M...
                + out.dY_R_dK_B/out.Y_R * out.W_Pi_R...
                + out.dY_C_dK_B/out.Y_C * out.W_Pi_C);
        
    % IC in eq (208)
    lambda_MR = out.dY_M_dK_B/out.Y_M-out.dY_R_dK_B/out.Y_R;
    lambda_CR = out.dY_C_dK_B/out.Y_C-out.dY_R_dK_B/out.Y_R;
    
    out.eq208 = -out.mu*(lambda_MR*out.I_M + lambda_CR*out.I_C) + ...
        out.mu*(  out.dY_M_dK_B/out.Y_M * out.I_Pi_M...
                + out.dY_R_dK_B/out.Y_R * out.I_Pi_R...
                + out.dY_C_dK_B/out.Y_C * out.I_Pi_C);
    
    % LC in eq (230)
    
    sum_Ci_M = out.xi_M*out.C_MM + out.xi_R*out.C_RM + out.xi_C*out.C_CM; 
    sum_Ci_C = out.xi_M*out.C_MC + out.xi_R*out.C_RC + out.xi_C*out.C_CC; 

    out.eq230 = -out.mu * (lambda_MR*sum_Ci_M + lambda_CR*sum_Ci_C)...
        + out.mu/out.K_B*out.C_Pi_B + out.mu/out.K_B*out.S_Pi_B;
    
    %% foc's for adding
    % base foc for V on (169)
    
    out.foc_V_row1 = sum(welfare_term_OB(gpops_fun,out).*...
                         out.V_prime./out.pi_w.*out.dH_dK_B.*out.dw);
    
    out.foc_V_row2 = -out.mu * sum(out.V_prime ./ out.h .* out.eta_prime ...
                         .* out.dH_dK_B .* out.dw);
    
    out.foc_V_row3 = out.mu * sum(out.l .* out.w .* out.pi_prime_w ...
                                  .* out.dH_dK_B .*...
                                  (  out.xi_M/out.Y_M * out.f_M./out.h ...
                                   + out.xi_R/out.Y_R * out.f_R./out.h ...
                                   + out.xi_C/out.Y_C * out.f_C./out.h) ...
                                  .* out.dw);
    
    out.foc_V_row4 = out.mu * sum((out.l.*out.w - out.c + out.b).* ...
                                  out.pi_prime_w./out.pi_w.*out.dH_dK_B.*out.dw);
    
    % row simplifies with quasi-linear utility
    out.foc_V_row5 = -out.mu * sum(out.V_prime .* out.dH_dK_B .* out.dw);
    
    out.foc_V = out.foc_V_row1 + out.foc_V_row2 + out.foc_V_row3 + ...
        out.foc_V_row4 + out.foc_V_row5;
    
    % base foc for l on (172)
    
    out.foc_l_row1 = out.mu * sum( out.l_prime .* out.w .* out.dH_dK_B ...
                                   .* (   out.xi_M/out.Y_M .*out.h_M./out.h ...
                                        + out.xi_R/out.Y_R .*out.h_R./out.h ...
                                        + out.xi_C/out.Y_C .*out.h_C./out.h)...
                                   .* out.dw);
    
    out.dmrs_dl = (1./epsilon)*1/xi*out.l.^(1./epsilon-1);
    
    out.foc_l_row2 = out.mu .* sum(out.l_prime./out.h .* ...
                                   (  out.w.*out.h.*(1 - out.mrs./out.w) ...
                                    - out.eta.* (  out.dmrs_dl.*out.l./out.w ...
                                                 + out.mrs./out.w ) )...
                                   .* out.dH_dK_B .* out.dw);
    
    out.foc_l = out.foc_l_row1 + out.foc_l_row2;
    
    %% bound check info
    out.check_lb_L_M = check_lb(1);
    out.check_lb_L_R = check_lb(2);
    out.check_lb_L_C = check_lb(3);
    out.check_lb_K_B_M = check_lb(4);
    out.check_lb_K_B_R = check_lb(5);
    out.check_lb_K_B_C = check_lb(6);
    out.check_lb_K_E = check_lb(7);
    out.check_lb_K_S = check_lb(8);
    out.check_lb_b = check_lb(9); 

    out.check_ub_L_M = check_ub(1);
    out.check_ub_L_R = check_ub(2);
    out.check_ub_L_C = check_ub(3);
    out.check_ub_K_B_M = check_ub(4);
    out.check_ub_K_B_R = check_ub(5);
    out.check_ub_K_B_C = check_ub(6);
    out.check_ub_K_E = check_ub(7);
    out.check_ub_K_S = check_ub(8);
    out.check_ub_b = check_ub(9);    
    
    
    %% analyze ranges
    tab.min_w_transf = output.result.solution.phase.time(1);
    tab.max_w_transf = output.result.solution.phase.time(end);
    tab.f_1 = out.f(1);
    tab.f_end = out.f(end);
    tab.f_M_1 = out.f_M(1);
    tab.f_M_end = out.f_M(end);
    tab.f_R_1 = out.f_R(1);
    tab.f_R_end = out.f_R(end);
    tab.f_C_1 = out.f_C(1);
    tab.f_C_end = out.f_C(end);
    tab.min_l_scaled = min(output.result.interpsolution.phase.control(:,1));
    tab.max_l_scaled = max(output.result.interpsolution.phase.control(:,1));
    tab.min_V_scaled = min(output.result.interpsolution.phase.state);
    tab.max_V_scaled = max(output.result.interpsolution.phase.state);
    tab.fval_scaled = output.result.objective;
    tab.L_M_scaled = output.result.solution.parameter(1,1);
    tab.L_R_scaled = output.result.solution.parameter(1,2);
    tab.L_C_scaled = output.result.solution.parameter(1,3);
    tab.K_B_M_scaled = output.result.solution.parameter(1,4);
    tab.K_B_R_scaled = output.result.solution.parameter(1,5);
    tab.K_B_C_scaled = output.result.solution.parameter(1,6);
    tab.K_E_scaled = output.result.solution.parameter(1,7);
    tab.K_S_scaled = output.result.solution.parameter(1,8);
    tab.b_scaled = output.result.solution.parameter(1,9);
    
    struct2table(tab)

    
    %% plot results
    figure(1)
    plot(out.w,out.mtax)
    title('Marginal Tax Rate')
    ylim([-0.3,1])
    xlabel('w')
    
    figure(2)
    plot(out.wtilde,out.mtax)
    title('Marginal Tax Rate')
    ylim([-0.3,1])
    xlabel('wtilde')
    
    figure(3)
    plot(out.y,out.mtax)
    title('Marginal Tax Rate')
    ylim([-0.5,1])
    xlabel('y')
    
    figure(4)
    plot(out.y,out.tax)
    title('Tax burden')
    xlabel('y')
    
    figure(5)
    plot(out.w,out.y)
    title('Income against w')
    xlabel('w')
    
    figure(6)
    plot(out.y,out.avg_tax)
    title('Average tax rate')
    xlabel('y')

    figure(8)
    plot(out.wtilde_noninterp,out.l_noninterp,'-x')
    title('l')
    xlabel('wtilde')    

    figure(9)
    plot(out.wtilde_noninterp,out.u_noninterp,'-x')
    title('V')
    xlabel('wtilde')        

    figure(10)
    plot(out.w_noninterp,out.l_noninterp,'-x')
    title('l')
    xlabel('w')    

    figure(11)
    plot(out.w_noninterp,out.u_noninterp,'-x')
    title('V')
    xlabel('w')         
    
    figure(12)
    plot(out.wtilde,out.f,...
         out.wtilde,out.f_M,...
         out.wtilde,out.f_R,...
         out.wtilde,out.f_C)
    legend('f','f_M','f_R','f_C')
    xlabel('wtilde')           
    
    figure(13)
    plot(out.wtilde,out.h,...
         out.wtilde,out.h_M,...
         out.wtilde,out.h_R,...
         out.wtilde,out.h_C)
    legend('h','h_M','h_R','h_C')
    xlabel('wtilde')     
    
    
    %% save results
    vars = rmfield(out,{'wtilde','w','u','l','c','mtax','y','tax','avg_tax',...
                        'wtilde_noninterp','w_noninterp','u_noninterp','l_noninterp',...
                        'l_orig','u_orig','eta_orig','wtilde_orig','w_orig',...
                        'f_M','f_R','f_C','f', ...
                        'h_M','h_R','h_C','h','F_M','F_R','F_C',...
                        'F','H_M','H_R','H_C','H',...
                        'h_M_cond_part','h_R_cond_part','h_C_cond_part',...
                        'h_cond_part','H_M_cond_part','H_R_cond_part','H_C_cond_part','H_cond_part',...
                       'eta','V_prime','dw','dh_M_h_dw','dh_R_h_dw','dh_C_h_dw','l_prime','pi_w','pi_prime_w',...
                       'c_prime','Pi_M','Pi_R','Pi_C','dPi_M_wh_dw','dPi_R_wh_dw','dPi_C_wh_dw',...
                       'dF_M_dK_B','dF_R_dK_B','dF_C_dK_B','dF_dK_B','dH_M_dK_B','dH_R_dK_B','dH_C_dK_B','dH_dK_B','mrs','phi_tilde_w','eta_prime','dMRS_l_dw','dmrs_dl'});

    vectors = keepfield(out,{'wtilde','w','u','l','c','mtax','y','tax','avg_tax',...
                        'f_M','f_R','f_C','f', ...
                        'h_M','h_R','h_C','h','F_M','F_R','F_C',...
                        'F','H_M','H_R','H_C','H',...
                        'h_M_cond_part','h_R_cond_part','h_C_cond_part',...
                        'h_cond_part','H_M_cond_part','H_R_cond_part','H_C_cond_part','H_cond_part',...
                       'eta','V_prime','dw','dh_M_h_dw','dh_R_h_dw','dh_C_h_dw','l_prime','pi_w','pi_prime_w',...
                       'c_prime','Pi_M','Pi_R','Pi_C','dPi_M_wh_dw','dPi_R_wh_dw','dPi_C_wh_dw',...
                       'dF_M_dK_B','dF_R_dK_B','dF_C_dK_B','dF_dK_B','dH_M_dK_B','dH_R_dK_B','dH_C_dK_B','dH_dK_B','mrs','phi_tilde_w','eta_prime','dMRS_l_dw','dmrs_dl'});
    
    
    vectors.q_B = vars.q_B .* ones(size(vectors.w));
    vectors.q_E = vars.q_E .* ones(size(vectors.w));
    vectors.q_S = vars.q_S .* ones(size(vectors.w));
    vectors.drop = vars.drop .* ones(size(vectors.w));
    
    dir_out = '../../../data/generated/matlab/OptimTax/'
toc    
end

function out = welfare_term_foc(self, var)
% numerically integrate over participation costs
% integration is done of the transformed variable chi(phi),
% which transforms the domain of integration from
% [self.phi_lb, var.phi_tilde_w] to [-1,1] and transforms
% phi to chi accordingly. In contrast to welfare in GPOPSFunctions,
% no multiplication with var.wscale in the last row, since here I
% will integrate over w's, not over log(w)
    
    n_nodes = length(self.nodes);

    chi = @(phi,lb,ub) ...
          tools.ChangeOfVariables.change_variable_minone_one(phi,lb,ub);
    
    dchi_dphi = @(phi,lb,ub) ...
        tools.ChangeOfVariables.change_variable_minone_one_deriv(phi,lb,ub);            

    fun = @(u,phi,lb,ub) Psi_prime(self,u - chi(phi,lb,ub)).*self.part_dist.g_M(chi(phi,lb,ub)).* ...
          dchi_dphi(phi,lb,ub);
    
    nodes = self.nodes(:)';

    phi_tilde = var.phi_tilde_w;
    phi_tilde(phi_tilde < self.phi_lb) = self.phi_lb;
    
    u_mat = repmat(var.u(:),1,n_nodes);
    lb_mat = self.phi_lb.*ones(size(u_mat));
    ub_mat = repmat(phi_tilde(:),1,n_nodes);
    nodes_mat = repmat(nodes,length(phi_tilde(:)),1);

    fun_at_nodes = fun(u_mat,nodes_mat,lb_mat,ub_mat);

    out = ((fun_at_nodes * self.weights(:)).* ...
           var.V_prime(:).*var.dF_dK_B(:));
end

function out = welfare_term_OB(self, var)
% numerically integrate over participation costs
% integration is done of the transformed variable chi(phi),
% which transforms the domain of integration from
% [self.phi_lb, var.phi_tilde_w] to [-1,1] and transforms
% phi to chi accordingly. In contrast to welfare in GPOPSFunctions,
% no multiplication with var.wscale in the last row, since here I
% will integrate over w's, not over log(w)
    
    n_nodes = length(self.nodes);

    chi = @(phi,lb,ub) ...
          tools.ChangeOfVariables.change_variable_minone_one(phi,lb,ub);
    
    dchi_dphi = @(phi,lb,ub) ...
        tools.ChangeOfVariables.change_variable_minone_one_deriv(phi,lb,ub);            

    fun = @(u,phi,lb,ub) Psi_prime(self,u - chi(phi,lb,ub)).*self.part_dist.g_M(chi(phi,lb,ub)).* ...
          dchi_dphi(phi,lb,ub);
    
    nodes = self.nodes(:)';

    phi_tilde = var.phi_tilde_w;
    phi_tilde(phi_tilde < self.phi_lb) = self.phi_lb;
    
    u_mat = repmat(var.u(:),1,n_nodes);
    lb_mat = self.phi_lb.*ones(size(u_mat));
    ub_mat = repmat(phi_tilde(:),1,n_nodes);
    nodes_mat = repmat(nodes,length(phi_tilde(:)),1);

    fun_at_nodes = fun(u_mat,nodes_mat,lb_mat,ub_mat);

    out = fun_at_nodes * self.weights(:);

end

function out = welfare_term_B_1(self, var)
% numerically integrate over participation costs
% integration is done of the transformed variable chi(phi),
% which transforms the domain of integration from
% [self.phi_lb, var.phi_tilde_w] to [-1,1] and transforms
% phi to chi accordingly. In contrast to welfare in GPOPSFunctions,
% no multiplication with var.wscale in the last row, since here I
% will integrate over w's, not over log(w)
    
    n_nodes = length(self.nodes);

    chi = @(phi,lb,ub) ...
          tools.ChangeOfVariables.change_variable_minone_one(phi,lb,ub);
    
    dchi_dphi = @(phi,lb,ub) ...
        tools.ChangeOfVariables.change_variable_minone_one_deriv(phi,lb,ub);            

    fun = @(u,phi,lb,ub) Psi(self,u - chi(phi,lb,ub)).*self.part_dist.g_M(chi(phi,lb,ub)).* ...
          dchi_dphi(phi,lb,ub);
    
    nodes = self.nodes(:)';

    phi_tilde = var.phi_tilde_w;
    phi_tilde(phi_tilde < self.phi_lb) = self.phi_lb;
    
    u_mat = repmat(var.u(:),1,n_nodes);
    lb_mat = self.phi_lb.*ones(size(u_mat));
    ub_mat = repmat(phi_tilde(:),1,n_nodes);
    nodes_mat = repmat(nodes,length(phi_tilde(:)),1);

    fun_at_nodes = fun(u_mat,nodes_mat,lb_mat,ub_mat);

    out = ((fun_at_nodes * self.weights(:)).* ...
           var.f(:));
end

function out = welfare_term_B_2(self, var)
% numerically integrate over participation costs
% integration is done of the transformed variable chi(phi),
% which transforms the domain of integration from
% [self.phi_lb, var.phi_tilde_w] to [-1,1] and transforms
% phi to chi accordingly. In contrast to welfare in GPOPSFunctions,
% no multiplication with var.wscale in the last row, since here I
% will integrate over w's, not over log(w)
    
    n_nodes = length(self.nodes);

    chi = @(phi,lb,ub) ...
          tools.ChangeOfVariables.change_variable_minone_one(phi,lb,ub);
    
    dchi_dphi = @(phi,lb,ub) ...
        tools.ChangeOfVariables.change_variable_minone_one_deriv(phi,lb,ub);            

    fun = @(u,phi,lb,ub) Psi(self,u - chi(phi,lb,ub)).*self.part_dist.g_M(chi(phi,lb,ub)).* ...
          dchi_dphi(phi,lb,ub);
    
    nodes = self.nodes(:)';

    phi_tilde = var.phi_tilde_w;
    phi_tilde(phi_tilde < self.phi_lb) = self.phi_lb;
    
    u_mat = repmat(var.u(:),1,n_nodes);
    lb_mat = self.phi_lb.*ones(size(u_mat));
    ub_mat = repmat(phi_tilde(:),1,n_nodes);
    nodes_mat = repmat(nodes,length(phi_tilde(:)),1);

    fun_at_nodes = fun(u_mat,nodes_mat,lb_mat,ub_mat);

    out = ((fun_at_nodes * self.weights(:)).* ...
           var.dF_dK_B(:));
end

function out = welfare_shift_M(self, var)
% numerically integrate over participation costs
% integration is done of the transformed variable chi(phi),
% which transforms the domain of integration from
% [self.phi_lb, var.phi_tilde_w] to [-1,1] and transforms
% phi to chi accordingly. In contrast to welfare in GPOPSFunctions,
% no multiplication with var.wscale in the last row, since here I
% will integrate over w's, not over log(w)
    
    n_nodes = length(self.nodes);

    chi = @(phi,lb,ub) ...
          tools.ChangeOfVariables.change_variable_minone_one(phi,lb,ub);
    
    dchi_dphi = @(phi,lb,ub) ...
        tools.ChangeOfVariables.change_variable_minone_one_deriv(phi,lb,ub);            

    fun = @(u,phi,lb,ub) Psi_prime(self,u - chi(phi,lb,ub)).*self.part_dist.g_M(chi(phi,lb,ub)).* ...
          dchi_dphi(phi,lb,ub);
    
    nodes = self.nodes(:)';

    phi_tilde = var.phi_tilde_w;
    phi_tilde(phi_tilde < self.phi_lb) = self.phi_lb;
    
    u_mat = repmat(var.u(:),1,n_nodes);
    lb_mat = self.phi_lb.*ones(size(u_mat));
    ub_mat = repmat(phi_tilde(:),1,n_nodes);
    nodes_mat = repmat(nodes,length(phi_tilde(:)),1);

    fun_at_nodes = fun(u_mat,nodes_mat,lb_mat,ub_mat);

    out = ((fun_at_nodes * self.weights(:)).* ...
           var.V_prime(:)./var.pi_w(:) .* var.Pi_M(:));
end

function out = welfare_shift_R(self, var)
% numerically integrate over participation costs
% integration is done of the transformed variable chi(phi),
% which transforms the domain of integration from
% [self.phi_lb, var.phi_tilde_w] to [-1,1] and transforms
% phi to chi accordingly. In contrast to welfare in GPOPSFunctions,
% no multiplication with var.wscale in the last row, since here I
% will integrate over w's, not over log(w)
    
    n_nodes = length(self.nodes);

    chi = @(phi,lb,ub) ...
          tools.ChangeOfVariables.change_variable_minone_one(phi,lb,ub);
    
    dchi_dphi = @(phi,lb,ub) ...
        tools.ChangeOfVariables.change_variable_minone_one_deriv(phi,lb,ub);            

    fun = @(u,phi,lb,ub) Psi_prime(self,u - chi(phi,lb,ub)).*self.part_dist.g_R(chi(phi,lb,ub)).* ...
          dchi_dphi(phi,lb,ub);
    
    nodes = self.nodes(:)';

    phi_tilde = var.phi_tilde_w;
    phi_tilde(phi_tilde < self.phi_lb) = self.phi_lb;
    
    u_mat = repmat(var.u(:),1,n_nodes);
    lb_mat = self.phi_lb.*ones(size(u_mat));
    ub_mat = repmat(phi_tilde(:),1,n_nodes);
    nodes_mat = repmat(nodes,length(phi_tilde(:)),1);

    fun_at_nodes = fun(u_mat,nodes_mat,lb_mat,ub_mat);

    out = ((fun_at_nodes * self.weights(:)).* ...
           var.V_prime(:)./var.pi_w(:) .* var.Pi_R(:));
end

function out = welfare_shift_C(self, var)
% numerically integrate over participation costs
% integration is done of the transformed variable chi(phi),
% which transforms the domain of integration from
% [self.phi_lb, var.phi_tilde_w] to [-1,1] and transforms
% phi to chi accordingly. In contrast to welfare in GPOPSFunctions,
% no multiplication with var.wscale in the last row, since here I
% will integrate over w's, not over log(w)
    
    n_nodes = length(self.nodes);

    chi = @(phi,lb,ub) ...
          tools.ChangeOfVariables.change_variable_minone_one(phi,lb,ub);
    
    dchi_dphi = @(phi,lb,ub) ...
        tools.ChangeOfVariables.change_variable_minone_one_deriv(phi,lb,ub);            

    fun = @(u,phi,lb,ub) Psi_prime(self,u - chi(phi,lb,ub)).*self.part_dist.g_C(chi(phi,lb,ub)).* ...
          dchi_dphi(phi,lb,ub);
    
    nodes = self.nodes(:)';

    phi_tilde = var.phi_tilde_w;
    phi_tilde(phi_tilde < self.phi_lb) = self.phi_lb;
    
    u_mat = repmat(var.u(:),1,n_nodes);
    lb_mat = self.phi_lb.*ones(size(u_mat));
    ub_mat = repmat(phi_tilde(:),1,n_nodes);
    nodes_mat = repmat(nodes,length(phi_tilde(:)),1);

    fun_at_nodes = fun(u_mat,nodes_mat,lb_mat,ub_mat);

    out = ((fun_at_nodes * self.weights(:)).* ...
           var.V_prime(:)./var.pi_w(:) .* var.Pi_C(:));
end

function out = Psi_prime(self, V)
    out = V.^(-self.r);
end

function out = Psi(self, V)
    out = (V.^(1-self.r)-1)./(1-self.r);
end